Métodos Numéricos 2018
Lista 01

Nesta lista, vamos trabalhar com o seguinte processo estocástico AR(1):

\[ z_t = \rho z_{t-1} + \varepsilon_t , \]

com \(\varepsilon_t \sim N(\mu,\sigma^2)\). Inicialmente, assumiremos que \(\mu=0\), \(\rho = 0.95\) e \(\sigma = 0.007\), calibração de Cooley & Prescott (1995).

Questão 1

Discretize o processo acima usando o método de Tauchen (1986). Use 9 pontos.

Ao discretizar um processo que assume valores contínuos, é necessário determinar o intervalo do grid, o número de pontos, e sua localização. Além disso, é necessário determinar probabiliades condicionais de transição entre os estados de forma consistente com o processo original.

Pelo método de Tauchen, o intervalo de um grid de \(N\) pontos é determinado por \([\mu - \theta_N, \mu + \theta_N]\), em que \(\theta_N = m \frac{\sigma}{\sqrt{1-\rho^2}}\) é o desvio padrão do processo estocástico, multiplicado por um parâmetro de escala \(m\). Os demais \(N-2\) pontos são distribuídos uniformemente no intervalo.

# Parâmetros Iniciais
N     = 9
rho   = 0.95
sigma = 0.007
mu    = 0
m     = 3
tauchen_grid = function(N, m, sigma, rho) {
  thetaN = m * sigma / sqrt(1-rho^2)  
  seq(-thetaN, thetaN, length.out = N)
}
thetaT = tauchen_grid(N, m, sigma, rho)
as.matrix(thetaT)

             [,1]
 [1,] -0.06725382
 [2,] -0.05044037
 [3,] -0.03362691
 [4,] -0.01681346
 [5,]  0.00000000
 [6,]  0.01681346
 [7,]  0.03362691
 [8,]  0.05044037
 [9,]  0.06725382

Seja \(\theta_{i,t}\) um estado no período \(t\), e \(\theta_{j, t+1}\) um estado no período \(t+1\). No método de Tauchen, as probabilidades de transição dos pontos interiores do grid são dadas por

\[ p_{i,j} = \Phi \left( \frac{\theta_j + \Delta \theta/2 - (1-\rho) \mu - \rho \theta_i}{\sigma} \right) - \Phi \left(\frac{\theta_j - \Delta \theta/2 - (1-\rho) \mu - \rho \theta_i}{\sigma}\right)\]

em que \(\Phi\) é a CDF da \(N(\mu, \sigma^2)\). Nos cantos do grid, temos que

\[ p_{i,1} = \Phi \left(\frac{\theta_1 - (1-\rho) \mu - \rho \theta_i + \Delta \theta/2}{\sigma} \right) \]

\[ p_{i,N} = 1 - \Phi \left(\frac{\theta_N - (1-\rho) \mu - \rho \theta_i - \Delta \theta/2}{\sigma} \right).\]

tauchen_P = function(grid, rho, sigma, mu) {
  N = length(grid)
  delta = (max(grid) - min(grid)) / (N-1)
  PT = matrix(NA_real_, N, N)
  for(i in 1:N) {
    for(j in 1:N) {
      if(j == 1) {
        PT[i,j] = pnorm( (grid[1] - (1-rho)*mu - rho*grid[i] + delta/2) / sigma )
      } else if(j==N) {
        PT[i,j] = 1 - pnorm( (grid[N] - (1-rho)*mu - rho*grid[i] - delta/2) / sigma )
      } else {
        PT[i,j] = 
          pnorm((grid[j] + delta/2 - (1-rho)*mu - rho*grid[i]) / sigma) -
          pnorm((grid[j] - delta/2 - (1-rho)*mu - rho*grid[i]) / sigma)
      }
    }
  }
  return(PT)
}
PT = tauchen_P(thetaT, rho, sigma, mu)
round(PT, 4)

        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]
 [1,] 0.7644 0.2347 0.0009 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
 [2,] 0.0592 0.7405 0.1997 0.0006 0.0000 0.0000 0.0000 0.0000 0.0000
 [3,] 0.0001 0.0747 0.7569 0.1679 0.0004 0.0000 0.0000 0.0000 0.0000
 [4,] 0.0000 0.0001 0.0931 0.7669 0.1396 0.0002 0.0000 0.0000 0.0000
 [5,] 0.0000 0.0000 0.0002 0.1147 0.7702 0.1147 0.0002 0.0000 0.0000
 [6,] 0.0000 0.0000 0.0000 0.0002 0.1396 0.7669 0.0931 0.0001 0.0000
 [7,] 0.0000 0.0000 0.0000 0.0000 0.0004 0.1679 0.7569 0.0747 0.0001
 [8,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0006 0.1997 0.7405 0.0592
 [9,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0009 0.2347 0.7644

rowSums(PT) # Linhas da matriz devem somar 1.

[1] 1 1 1 1 1 1 1 1 1

Questão 2

Discretize o processo acima usando o método de Rouwenhorst. Use 9 pontos.

A definição do grid no método de Rouwenhorst é um caso particular do método de Tauchen, com \(m = \sqrt{N-1}\).

thetaR = tauchen_grid(N, m = sqrt(N-1), sigma, rho)
thetaR

[1] -0.06340751 -0.04755564 -0.03170376 -0.01585188  0.00000000  0.01585188  0.03170376  0.04755564  0.06340751

A matriz de transição \(P\), porém, é obtida recursivamente.

rouwen_P = function(N, rho) {
  p = (1+rho)/2
  P_init = matrix(c(p, 1-p, 1-p, p), 2, 2, byrow = TRUE)
  for(i in 3:N) {
    PR = 
      rbind(cbind(p*P_init, 0), 0) +
      rbind(cbind(0, (1-p)*P_init), 0) +
      rbind(0, cbind((1-p)*P_init, 0)) +
      rbind(0, cbind(0, p*P_init))
    PR = PR / rowSums(PR) # normalização
    P_init = PR
  }
  
  return(PR)
}
PR = rouwen_P(N, rho)
round(PR, 4)

        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]
 [1,] 0.8167 0.1675 0.0150 0.0008 0.0000 0.0000 0.0000 0.0000 0.0000
 [2,] 0.0209 0.8204 0.1469 0.0113 0.0005 0.0000 0.0000 0.0000 0.0000
 [3,] 0.0005 0.0420 0.8231 0.1261 0.0081 0.0003 0.0000 0.0000 0.0000
 [4,] 0.0000 0.0016 0.0630 0.8247 0.1051 0.0054 0.0001 0.0000 0.0000
 [5,] 0.0000 0.0001 0.0032 0.0841 0.8253 0.0841 0.0032 0.0001 0.0000
 [6,] 0.0000 0.0000 0.0001 0.0054 0.1051 0.8247 0.0630 0.0016 0.0000
 [7,] 0.0000 0.0000 0.0000 0.0003 0.0081 0.1261 0.8231 0.0420 0.0005
 [8,] 0.0000 0.0000 0.0000 0.0000 0.0005 0.0113 0.1469 0.8204 0.0209
 [9,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0008 0.0150 0.1675 0.8167

Questão 03

Simule o processo contínuo para 10000 períodos. Faça o mesmo para os processos discretizados (lembre-se de usar as mesmas realizações para os choques). Compare os caminhos para cada processo (gráficos serão úteis aqui). Se eles não estiverem muito próximos, utilize mais pontos.

set.seed(96452)
time_span = 10000
eps = rnorm(time_span, mu, sigma)
# Processo contínuo
zc = rep(0, time_span)
for(t in 2:time_span)
  zc[t] = rho * zc[t-1] + eps[t]
# Funcão para gerar processos
gen_process = function(grid, P, eps, mu, sigma) {
  time_span = length(eps)
  z = rep(0, time_span)
  for(t in 2:time_span) {
    i = which(grid == z[t-1])
    j = sum(cumsum(P[i, ]) < pnorm(eps[t], mu, sigma)) + 1
    z[t] = grid[j]
  }
  return(z)  
}
zt = gen_process(thetaT, PT, eps, mu, sigma)
zr = gen_process(thetaR, PR, eps, mu, sigma)

Os gráficos a seguir exibem o processo contínuo e as discretizações de Tauchen e de Rouwenhorst para 10.000 períodos e para os 500 primeiros períodos:

Visualmente, a discretização de Tauchen funciona melhor com os parâmetros dados. A discretização de Tauchen captura melhor a amplitude do processo sem perder precisão. Ambos apresentam, contudo um RMSE semelhante em relação ao processo contínuo:

RMSE = function(x,y) sqrt(mean((x-y)^2))
cat("RMSE:\n")
cbind("Tauchen" = RMSE(zc,zt), "Rouwenhorst" = RMSE(zc,zr))

Questão 4

Estime processos AR(1) com base nos dados simulados, tanto a partir do Tauchen quanto o de Rouwenhorst. Quão próximo eles estão do processo gerador de dados real? Se eles não estiverem muito próximos, utilize mais pontos.

fit_zt = lm(zt ~ shift(zt) + 0)
fit_zr = lm(zr ~ shift(zr) + 0)
coefs = c(fit_zt$coefficients,
          fit_zr$coefficients)
sigmas = c(summary(fit_zt)$sigma,
           summary(fit_zr)$sigma)
tbl = data.frame(method = c("Tauchen 9 pts", 
                            "Rouwenhorst 9 pts"), 
                 rho = coefs, sigma = sigmas)
tbl

Pelo método de Tauchen com 9 pontos no grid, estimamos \(\hat{\sigma} = 0.008\), valor relativamente distante do \(\sigma\) verdadeiro. Aumentando o número de pontos no grid para 51, obtemos uma estimativa semelhante à obtida para o método de Rouwenhorst com 9 pontos.

thetaT = tauchen_grid(N=51, m=3, sigma, rho)
PT = tauchen_P(thetaT, rho, sigma, mu)
ztl = gen_process(thetaT, PT, eps, mu, sigma)
fit_ztl = lm(ztl ~ shift(ztl) + 0)
rbind(tbl, cbind(method = "Tauchen 25 pts", rho = fit_ztl$coefficients, sigma = summary(fit_ztl)$sigma))

Questão 5

Vamos refazer os exercícios acima com \(\rho = 0.99\).

O método de Rouwenhorst funciona melhor para \(\rho = 0.99\) com 9 pontos no grid. Obtemos os seguintes valores de \(\hat{\rho}\) e \(\hat{\sigma}\):